---
title: "replication file"
author: "Thalia Gerzso & Nicolas van de Walle
date: "11/22/2021"
output: html_document
---

```{r}
#Required packages to run the survival analysis 
library(ggplot2)
library (dplyr)
library (tidyverse)
library(haven)
library(margins)
library(magrittr)
library(survival)
library(survminer)
library(GGally)
```

```{r}
#load the dataset
df_all <- read_csv("dflegi.csv")
```



```{r}
# Survival analysis Regression
res.cox14 <- coxph(formula = Surv(OBS, CENS) ~  TERMLIMIT + DEFEAT0 + INCPOP  +  GDP + NATRESOURCES + CAM1 + EFFPARTY + PARL + AVRFH + PR + DENSITY + APPEC90 + CONST, data = df_all)
summary(res.cox14)
```

```{r}
#Appendix 1 Shared Frailty with Gamma Fit Distribution
frail1 <- coxph(Surv(OBS, CENS) ~ TERMLIMIT  +  DEFEAT0 + INCPOP + GDP + NATRESOURCES + CAM1 + EFFPARTY + PARL + AVRFH + PR + DENSITY +  APPEC90 + CONST + frailty.gamma(COUNTRY), df_all)
summary(frail1)


#Appendix 2: Schoenfeld residual tests 
test.ph <- cox.zph(res.cox14)
test.ph
res.cox16 <- coxph(formula = Surv(start, stop,  CENS) ~ TERMLIMIT  + DEFEAT0 + INCPOP + GDP + NATRESOURCES + CAM1 + EFFPARTY + PARL + INCPOP:stop + AVRFH + PR + DENSITY + APPEC90 + CONST , data = df_all)
test.ph2 <- cox.zph(res.cox16)
test.ph2

#Appendix 3: regression removing democratic states
df_dem <-
  filter (df_all, AVRFH > 5)

res.cox15 <- coxph(formula = Surv(OBS, CENS) ~ TERMLIMIT + DEFEAT0 + INCPOP + GDP + NATRESOURCES + CAM1 + EFFPARTY  + PARL + AVRFH + PR + DENSITY + APPEC90 + CONST, data = df_dem)
summary(res.cox15)

#Appendix 4: Correlation between executive power measures
df3 <- select(df_all, VDEMEXE, PRESPOW1, PRESPOW2, CPP, XCONSTAVR,TERMLIMIT)
df4 <-df3 %>% 
  select(VDEMEXE, PRESPOW1, PRESPOW2, CPP, XCONSTAVR,TERMLIMIT) %>% 
  rename ("VDEM" = VDEMEXE,
          "PP1" = PRESPOW1,
          "PP2" = PRESPOW2, 
          "CPI" = CPP, 
          "Polity" = XCONSTAVR,
          "TL" = TERMLIMIT
        )
ggcorr(df4, label = TRUE)
```
```{r}
# plot comparing sub-saharan Africa, Eastern Europe, and Latin America
df1 <- read_csv("World-legi.csv")


df2 <- aggregate(df1, 
                  by= list(df1$Year, df1$Region), 
                  FUN = function(x) mean(x, na.rm = TRUE))
ggplot(df2, aes(x = Year, y = Size,
                 color = factor(Region2))) +
  geom_line() +
   scale_colour_discrete(name  ="Region",
                           breaks=c("1", "3","4"),
                           labels=c("Latin America", "Eastern Europe", "Sub-Saharan Africa")) +
  theme_bw() + 
  labs(y="Legislature size", x= "Years" ) +
  ggsave ("region.pdf", width = 5.29, height = 2.51)
  
 ggplot(df2, aes(x = Year, y = Low,
                 color = factor(Region2))) +
  geom_line() +
   scale_colour_discrete(name  ="Region",
                           breaks=c("1", "3","4"),
                           labels=c("Latin America", "Eastern Europe", "Sub-Saharan Africa")) +
  theme_bw() + 
  labs(y="Legislature size", x= "Years" )

 # descriptive statistics comparing sub-saharan Africa, Eastern Europe, and Latin America 
 df1 %>% 
 filter( Region == "Sub-Saharan Africa", Change == 1) %>% 
  select (Size, Change, Bicam, Inc, Dec, Low ) %>% 
   summary (Region)

df1 %>% 
 filter( Region == "Sub-Saharan Africa") %>% 
  count (Change)

df1 %>% 
  filter (Region == "Sub-Saharan Africa", Bicam == 1) %>% 
  count (Country)

df1 %>% 
 filter( Region == "Sub-Saharan Africa",Change == 1) %>% 
  count(Country)

